Optimal system size for complex dynamics in random neural networks near criticality 
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In this Letter, we consider a model of dynamical agents coupled through a random connectivity 
matrix, as introduced in [Sompolinsky et. al, 1988] in the context of random neural networks. 
It is known that increasing the disorder parameter induces a phase transition leading to chaotic 
dynamics. We observe and investigate here a novel phenomenon in the subcritical regime : the 
probability of observing complex dynamics is maximal for an intermediate system size when the 
disorder is close enough to criticality. We give a more general explanation of this type of system 
size resonance in the framework of extreme values theory for eigenvalues of random matrices. 
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Most complex systems involve a large number of inter- 
acting elements. When a detailed knowledge of the con- 
nections' properties is lacking, statistical modeling of the 
connectivity is often introduced giving rise to the study of 
dynamical processes evolving on random networks. Typ- 
ical examples range from condensed matter physics pQ , to 
biology [HG3, and social sciences [I]. The hundreds of bil- 
lions of heterogeneous synaptic connections between neu- 
rons in the brain constitute a paradigmatic example of a 
disordered system of dynamical agents coupled through a 
random connectivity matrix. To understand the impact 
of this heterogeneity on brain dynamics, a random neural 
network model has been introduced in [2J: 

For 1 < i < n, ^ = -x* + ^ Jij</>(xj) (1) 

where J = (Jy) is a Ginibre [5] random matrix with 
Gaussian i.i.d entries such that E[J#] = and E[J?-] = 
a 2 /n, and where 4>{-) is a smooth odd sigmoid function 
with unit slope at the origin </>'(0) = 1. Notice that zero 
is a trivial equilibrium of (JlJ. 

It is known since [2] (see also [6]) that in the limit 
ri — > oo one can derive a mean field description of the 
above model and an analysis of the mean field equations 
[2] leads to the two following regimes : 

(i) if a < 1 the only equilibrium point is zero and 
attracts "all" trajectories ; 

(ii) if a > 1 there exist infinitely many equilibria and 
limit cycles but none of them is stable. The only stable 
attractors are chaotic in the sense of positive Lyapunov 
exponents. 

Therefore, in this system, the level of disorder induces 
a phase transition. However, as stressed in [2J, the pic- 
ture is different for finite size systems where, close to the 
mean field phase transition, several stable and chaotic 
attractors may co-exist. The stable attractors are either 
non zero equilibrium points or limit cycles. 

Based on a combination of experimental and theoreti- 
cal considerations, it has been argued in 0[S] that living 
neuronal networks, as well as many other complex sys- 



tems [S], may evolve in the vicinity of the phase transi- 
tion also called the edge of chaos [10 . In this region the 
extreme eigenvalues of J play a crucial role in the stabil- 
ity of finite systems, suggesting that the behavior of the 
system will not be properly described by the mean-field 
equations. 

Therefore, a fine knowledge of the statistical proper- 
ties of the eigenvalues distribution of random matrices is 
fundamental to characterize the behavior of various dis- 
ordered systems, in particular in terms of stability prop- 
erties [11] . One of the most important results concerning 
non-hermitian matrices is Girko's circular law |12j and 
has been recently proven universal in [13] . This theo- 
rem states that the empirical distribution of the com- 
plex eigenvalues for n x n random matrices with i.i.d 
centered coefficients of variance n~ l , converges to the 
uniform measure on the unit disk. In the case of finite 
size matrices, despite most eigenvalues being inside the 
unit circle, there is a nonzero probability that eigenvalues 
close to the border have a modulus slightly larger than 
1, leading to unstable modes from the stability point of 
view. 

Our aim is to study the probability of observing spon- 
taneous activity in the network near the phase tran- 
sition, specifically for values of disorder slightly below 
criticality. As n —¥ oo, this probability shall converge 
to zero, but as we will show, this convergence is not 
monotonous, this probability is actually reaching a max- 
imum for an intermediate value of n. We will first inves- 
tigate this question directly on system ([I]) numerically, 
and then provide a deeper theoretical explanation of this 
phenomenon of system size resonance through the study 
of extreme eigenvalues of random matrices. 

To determine wether the network is activated or not 
we compute the maximal Lyapunov exponent A of the 
attractors. Depending on the value of A one can distin- 
guish three cases: (i) if A < the attractor is a fixed 
point, the network is not activated, (ii) if A = the 
attractor is a limit cycle, the network is activated and 
exhibits regular activity, and (iii) if A > the attractor 
is chaotic, the network is activated and exhibits irregu- 
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lar activity. In terms of A, the probability of observing 
spontaneous activity, whether it is chaotic or periodic, is 
given by P[A > 0]. 

In the numerical simulations, an estimation of the max- 
imal Lyapunov exponent of the attractors A is computed 
using the variational equations as shown in [14] . Starting 
from a random initial condition, we evolve a reference so- 
lution (xi(t))i<i< n according to ([I]) and a normalized dif- 
ference vector on the tangent plane (6(x, i)i)i<i<n driven 
by the Jacobian of 0. At each time step the mod- 
ulus of S is normalized to a small quantity 5q so that 
(x(t)i + 5(x, t)i)i<i< n is always in the vicinity of the ref- 
erence solution. After some iterations the reference orbit 
will be close to one of the attractors of the system and the 
difference vector will point in the direction corresponding 
to the maximal Lyapunov exponent of the reference tra- 



jectory. Then, we estimate A = * J2ilT 1°S jj, where T 
is the sampling time and dt is the time step. To circum- 
vent approximations error in the case of periodic attrac- 
tors, we set A = when periodicity is detected directly 
on the trajectories [15] . 

Numerical results are shown in Fig. [I] The main unex- 
pected phenomenon is that the probability of observing 
spontaneous activity for a < 1, close enough to 1, dis- 
plays a maximum for an intermediate value of the net- 
work size n (Fig.[TJa.). Surprisingly, this probability first 
increases with n until an optimal system size n(<j), and 
then decreases to zero as expected from [2]. Moreover, 
as a gets closer to 1, both the optimal size n(a) and the 
probability of spontaneous activity at n = n(a) increase. 
In the limit case a = 1 we do not observe a maximum 
and P[A > 0] tends to 1 as n becomes larger. 

The estimations of the probabilities corresponding to 
limit cycles P[A = 0] (Fig. [l]b.) and chaotic oscillations 
P[A > 0] (Fig. [l]c.) have a similar behavior. The main 
difference is that the maxima for P[A > 0] are reached at 
larger values of n than the ones for P[A = 0] . As shown in 
Fig. [TJd., this difference means that the larger the system 
is, the more likely is that the spontaneous activity takes 
the form of chaotic attractors. 

From numerical simulations, we have shown the exis- 
tence of an intermediate system size which maximizes the 
probability of observing complex dynamics. This myste- 
rious phenomenon has not been described previously and 
we intend to provide a theoretical explanation through 
the lens of random matrix theory. Indeed, as mentioned 
in [2] , the phase transition can be inferred from the study 
of the linearized problem around the trivial equilibrium: 
in the limit n = oo, linear stability is lost when a > 1 
as a direct consequence of the circular law. Moreover, if 
(Afc)i</c<„ denote the complex eigenvalues of J, and if 
one defines 



pi = max Re(\k), 

Kk<n 



(2) 



then pi < 1 implies that the trivial equilibrium is globally 
attractive. Therefore, the natural quantity to study is the 



probability p n that the trivial equilibrium zero is linearly 
stable. Such a question has a direct translation in terms 
of random matrix theory: 



Pn =P[ P i< 1] 



(3) 



First, let us remark the fact that if er < 1 then p„ 1 
when Tt — y oo . We will show that this probability has 
a minimum for an intermediate matrix size, which in- 
creases to +oo when a — > 1~ . In terms of the original 
system (JTJ), this relates to the fact that the probability 
of observing a spontaneous oscillations or complex dy- 
namics when a is close to 1 _ displays a maximum for an 
intermediate system size. The probability p n is not ex- 
actly equal to P[A < 0], but one only has the inequality 
Pn < P[A < 0]. Indeed, it is possible to have p\ > 1 
and still A < due to the existence non trivial attracting 
equilibrium points. However, we expect both quantities 
to have similar properties. 

To study pi , we need to build upon the extreme value 
theory for the eigenvalues random matrices. Contrary to 
classical extreme value theory, there is a deep dependance 
structure among the eigenvalues, which is well described 
within the framework of determinental point processes. 
In fact, the theory of extremal eigenvalues of complex 
Gaussian matrices has been investigated recently in [16] . 

Based on these results, we explain the optimal size 
phenomenon in the case of complex matrices. In the case 
of real Gaussian matrices, the situation is more delicate 
and we only provide numerical results in Fig. [2] showing 
the same phenomenology. 

More precisely, for a > 0, we consider {A^(cr)} i<k< n 
the complex eigenvalues of a n x n random matrix J 
such that the coefficients Jjj are i.i.d complex Gaussian 
random variables on C = M. 2 with Jy « A/"(0, §— -^2)- 
We denote the maximum of the real parts of the eigen- 
values of J by: 



Kk<n 



(4) 



We are interested in the probability that all the eigenval- 
ues have a real part smaller than 1: 



p n (a) :=p[aW» 



< 1 



(5) 



and its complementary probability q n {o') = 1 — Pn(cr), 
that is the probability that there exists at least one eigen- 
value with a real part larger than 1. Our main theoretical 
result is that, for any a £ (0, 1), the integer n*(cr) such 
that q n *(a) — sup q n (o ) satisfies: 



lim n*(a) 

7 — y 1 



—(—00 



Moreover, q n (a) can be made arbitrary close to 1: 



lim q n .r a ) 
7— >i~ 



(G) 



(7) 
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FIG. 1. Numerical estimation of the probability of observing (a) spontaneous activity (b) limit cycles (c) chaotic trajectories and 
(d) chaotic attractor given spontaneous activity, as a function of n and for different values of a £ {0.95, 0.96, 0.97, 0.98, 0.99, 1}. 
For each matrix size and each value of a < 1, 4.10 4 realizations of the random matrix have been analyzed. Numerical estimations 
of the maximal Lyapunov exponent were performed after cutting a first period of transient dynamics. 



Heuristically, this property stems from a competition be- 
tween (i) the fact that when the number of eigenvalues n 
increases, the maximal real part tends to be higher, as in 
classical extreme value theory, and (ii) the convergence 
of the spectral density to the unit disk that implies a 
concentration of \mlx(cr) close to a < 1. 

We now provide a justification of these results, based 
on |16| , where the asymptotic law for large n of the largest 
real part of eigenvalues for a class non-Hcrmitian random 
matrix is studied, developing extreme value theory for 
determinental processes. Applying Theorem 2.5 of [16], 
we obtain: 



lim I 

n— >oo 



A (n) (1) - 1< 



l 



zt + C r< 



2Vnlogn 
for any t G K, with 

1 fW) ^HHn))-^/^) 



2^2 V n ^/nln(n) 
From this result, we claim that for any 5 € (0, 1): 



lim P 

n— >oo 



\%l(l)<l + (l -6)c n 







(8) 



(9) 



(10) 



Indeed, let 5 G (0,1) and e > fixed, and t, < such 
that e~ e ' < e. Then from (JsJ) there exists uq > 1 such 
that for all n > n : 



A^L(l)<l + Cnfl+ _ 1— ) 
V 2c n ^n\ognJ 



< e 



(11) 



Moreover, as c n \/nlogn — > 00 when n — > 00, one can 
choose tiq such that — S < t e /(2c n y/n log n). Then, by 
monoticity of the distribution function, we deduce: 

A£L(1)< 1 + ^(1-5) 

A&L(1)<1+Cnl I 



< 



2c n y / n logn 



< e 



which proves claim ( 10 1 



Now we want to show that p n (a) can be made arbitrarily 
small when a — > 1~. First, since the law of Xmax{<y) is 
the same as the law of a\ml x (l), we observe that: 

P [aW x (1)<1 + (1-5)c„ 

= F fe» < 1 + (7(1 - S)c n -(1-ct) := /„ 
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We know that f n — ¥ when n — > oo, so there exists 
m > 1 (independent of a) such that for all n > ni, 
/„ < e. Moreover, since c„ — > 0, we can choose h(a) 
such that for all n > n(cr), er(l — S)c n — (1 — cr) > 0, say 
we choose it such that (1 — 5)cf b *^ ~ 2(1 — cr)/cr. As 
n(cr) — » oo when a — > 1 — , there exist oo € (0, 1) such that 
"(o'o) > n o- Then, by monotonicity of the distribution 
function, wc deduce that: 



Pn( CTo )( CT o) < /n(a ) < 6 



(12) 



As a consequence, we deduce that n* (a) — > oo when cr — > 
1~ as stated in |6]) and that g„»( CT ) -> 1 as stated in Q. 

As far as the real Ginibre ensemble is concerned, the 
above theoretical analysis does not apply directly. In- 
deed, the spectral joint probability density in this case 
lacks the determinental structure, which is key in the 
work of |16j . However, very recently extreme value theory 
for the spectral radius of real matrices has been studied 
[T7] . suggesting that a similar result could be obtained 
in this context. In fact, numerical simulations confirm 
that the same phenomenon occurs. In Fig. [2j a numeri- 
cal estimation of q n {o) is computed and shows that for a 
close enough to 1, q n {<r) first increases then decreases as 
a function of n, similarly to the maximum observed on 
the maximal Lyapunov exponent. 

This phenomenon of system size resonance may appear 
in a wider class of problems, involving the maximum of 
a set n convergent random variables. It is instructive to 
discuss this idea on a toy example. Consider a sequence 

of families of real i.i.d random variables ( ) and 

V / i<j< n 

We assume that 



denote M n = max(A, 

(n) 



(«). 



1 < i < n) 



< x] = a n = 1 — b n is increasing to 1 when n — > oo 
for a given ieR. By the i.i.d assumption, one writes: 



»[m„ > x] = i - (i - b n y 



(13) 



We deduce that if b n = o(n _1 ) then P[M„ > x] -> when 
n — > oo. Now we can wonder whether this convergence 
is monotonic or if there exists an intermediate value of n 
for which V[M n > x] is maximal. This can be answered 
by differentiating ( 13 1 with respect to n (considered as a 
real variable) and looking n such that: 



a n log(a„) + na' n = 



(14) 



We find for instance that if a n = 1 — e~ Kn with k < 1, 
then P[Af„ > a;] is maximal for n around k" 1 . 

The real parts of the eigenvalues of J satisfy similar 
assumptions as the toy model, except the strongest one 
that is independence. This difficulty has been solved us- 
ing the approach of determinental processes in [16] , which 
was the starting point of our analysis. 

Moreover, one can view the results on the Lyapunov 
exponent as a consequence of this general principle. 
There is a competition between two opposite phenomena 
as the system size increases. On the one hand, consis- 
tent with an extreme value behavior, increasing system 
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FIG. 2. Real Ginibre Ensemble : Numerical estimation 
of q n (cr) as a function of n, for different values of a = 
0.97,0.98,0.99,1 (from bottom to top line). A total of 5.10 5 
Monte-Carlo runs were computed for each different matrix 
size and value of a. 



size increases the likelihood of obtaining a large maxi- 
mal Lyapunov exponent, since one looks at the maximum 
on n variables. On the other hand, as n increases, self- 
averaging principle drives the system towards the mean 
field behavior which converges to the trivial zero solution 
for a < 1. This competition results in the emergence of 
an intermediate system size which enhances the proba- 
bility of observing spontaneous activity. 

In terms of perspectives, our result may shed a new 
light on the behavior of modular networks, composed by 
clusters of randomly connected agents and sparse con- 
nections between clusters. Indeed, each cluster size may 
be viewed as a parameter that controls the propensity to 
generate complex dynamics within each cluster. 

Finally, although we have carried out the numerical 
investigation of the maximal Lyapunov exponent only 
on the random neural network model 0, the theoreti- 
cal analysis of Xmlx using random matrix theory is much 
more general and suggests that the system resonance phe- 
nomenon may be observed in various complex systems, 
and in particular in disordered systems close to the phase 
transition. 
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